####################################################################
# Replication code for
#   Hirose, Imai, and Lyall, "Can civilian attitudes predict insurgent violence?"
#
# Code to conduct in-sample descriptive & regression analysis 
#
# Author: Kentaro Hirose
####################################################################



########################################################################
library(endorse)
library(MASS)
library(mgcv)
library(sandwich)
library(robust)
library(Hmisc)
library(bayesm)
library(hglm)
library(lme4)
########################################################################
load("village.data.RData") ### in-sample village data
########################################################################
### in-sample violence (against ISAF)
load("pre.ied.attack.RData")
load("pre.ied.found.RData")
load("pre.nonied.RData")
load("post.ied.attack.RData")
load("post.ied.found.RData")
load("post.nonied.RData")
########################################################################
### in-sample violence (against civilians)
load("pre.ied.attack.civilian.RData")
load("pre.nonied.civilian.RData")
load("post.ied.attack.civilian.RData")
load("post.nonied.civilian.RData")
########################################################################




########################################################################
### Figure 2: 15 km, 150 days (effects of relative support)
########################################################################
k <- 15
t <- 5
pre.A <- pre.ied.attack[[t]][, k]
pre.B <- pre.ied.found[[t]][, k]
pre.C <- pre.nonied[[t]][, k]
post.A <- post.ied.attack[[t]][, k]
post.B <- post.ied.found[[t]][, k]
post.C <- post.nonied[[t]][, k]
########################################################################
J <- 100
mean1 <- mean2 <- mean3 <- mean4 <- mean5 <- mean6 <- mean7 <- mean8 <- mean9 <- mean10 <- mean11 <- mean12 <- rep(NA, J)
lower1 <- lower2 <- lower3 <- lower4 <- lower5 <- lower6 <- lower7 <- lower8 <- lower9 <- lower10 <- lower11 <- lower12 <- rep(NA, J)
upper1 <- upper2 <- upper3 <- upper4 <- upper5 <- upper6 <- upper7 <- upper8 <- upper9 <- upper10 <- upper11 <- upper12 <- rep(NA, J)
########################################################################
X<- village.data$support.isaf - village.data$support.taliban
aid <- village.data$aid
base <- village.data$base_3km
MIN <- min(X, na.rm = TRUE)
MAX <- max(X, na.rm = TRUE)
R <- seq(MIN, MAX, length.out = J)
########################################################################
### robust se ###
Y <- post.A
Y.lag <- pre.A
out <- lm(Y ~ X + Y.lag + base + aid)
out.poi <- glm(Y ~ X + Y.lag + base + aid, family = "poisson")
out.neg <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 100))
coef.out <- coef(out)
vcov.out <- vcovHC(out)
set.seed(1000)
beta <- mvrnorm(1000, mu = coef.out, Sigma = vcov.out)
mean.Y.lag <- mean(Y.lag, na.rm = TRUE)
mean.base <- mean(base, na.rm = TRUE)
median.aid <- median(aid, na.rm = TRUE)
Z <- cbind(1, R, mean.Y.lag, mean.base, median.aid)  
mu <- beta %*% t(Z)
for (j in 1:J){
  mu.j <- mu[, j]
  mean1[j] <- mean(mu.j)
  upper1[j] <- quantile(mu.j, prob = 0.975)
  lower1[j] <- quantile(mu.j, prob = 0.025)    
}

Y <- post.B 
Y.lag <- pre.B 
out <- lm(Y ~ X + Y.lag + base + aid)
out.poi <- glm(Y ~ X + Y.lag + base + aid, family = "poisson")
out.neg <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 100))
coef.out <- coef(out)
vcov.out <- vcovHC(out)
set.seed(1000)
beta <- mvrnorm(1000, mu = coef.out, Sigma = vcov.out)
mean.Y.lag <- mean(Y.lag, na.rm = TRUE)
mean.base <- mean(base, na.rm = TRUE)
median.aid <- median(aid, na.rm = TRUE)
Z <- cbind(1, R, mean.Y.lag, mean.base, median.aid)   
mu <- beta %*% t(Z)
for (j in 1:J){
  mu.j <- mu[, j]
  mean2[j] <- mean(mu.j)
  upper2[j] <- quantile(mu.j, prob = 0.975)
  lower2[j] <- quantile(mu.j, prob = 0.025)    
}

Y <- post.C 
Y.lag <- pre.C 
out <- lm(Y ~ X + Y.lag + base + aid)
out.poi <- glm(Y ~ X + Y.lag + base + aid, family = "poisson")
out.neg <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 100))
coef.out <- coef(out)
vcov.out <- vcovHC(out)
set.seed(1000)
beta <- mvrnorm(1000, mu = coef.out, Sigma = vcov.out)
mean.Y.lag <- mean(Y.lag, na.rm = TRUE)
mean.base <- mean(base, na.rm = TRUE)
median.aid <- median(aid, na.rm = TRUE)
Z <- cbind(1, R, mean.Y.lag, mean.base, median.aid)   
mu <- beta %*% t(Z)
for (j in 1:J){
  mu.j <- mu[, j]
  mean3[j] <- mean(mu.j)
  upper3[j] <- quantile(mu.j, prob = 0.975)
  lower3[j] <- quantile(mu.j, prob = 0.025)    
}

########################################################################
### computing change in attacks
########################################################################
RR <- c(-2.5, 0.5)
Z <- cbind(1, RR, mean.Y.lag, mean.base, median.aid)   
mu <- beta %*% t(Z)
diff <- mu[,2] - mu[,1]
quantile(diff, prob = 0.025)
mean(diff)
quantile(diff, prob = 0.975)
########################################################################
pdf(file = "effect_15km_5month.pdf", width = 12, height = 4)

par(mfrow = c(1,3))

#################################
X<- village.data$support.isaf - village.data$support.taliban
XLAB = "Relative level of ISAF support"      
YLAB <- "Number of insurgent attacks"
ADJ.FACTOR <- 50
#################################

MIN <- min(lower1, lower2, lower3)
MAX <- max(upper1, upper2, upper3)

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=2.2
cex.sub=x

LOWER <- lower1
MEAN <- mean1
UPPER <- upper1
MAIN <- "IED attacks"
ADJ <- abs(MAX - MIN)/ADJ.FACTOR
plot(R, MEAN, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN - ADJ, MAX), col = 0, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(R, LOWER, lty = 2, lwd = 1)
lines(R, MEAN, lty = 1, lwd = 3)
lines(R, UPPER, lty = 2, lwd = 1)
segments(x0 = X, x1 = X, y0 = MIN - 100, y1 = MIN, lwd = 0.2)

LOWER <- lower2
MEAN <- mean2
UPPER <- upper2
MAIN <- "IED found"
ADJ <- abs(MAX - MIN)/ADJ.FACTOR
plot(R, MEAN, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN - ADJ, MAX), col = 0, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(R, LOWER, lty = 2, lwd = 1)
lines(R, MEAN, lty = 1, lwd = 3)
lines(R, UPPER, lty = 2, lwd = 1)
segments(x0 = X, x1 = X, y0 = MIN - 100, y1 = MIN, lwd = 0.2)

LOWER <- lower3
MEAN <- mean3
UPPER <- upper3
MAIN <- "Non-IED attacks"
ADJ <- abs(MAX - MIN)/ADJ.FACTOR
plot(R, MEAN, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN - ADJ, MAX), col = 0, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(R, LOWER, lty = 2, lwd = 1)
lines(R, MEAN, lty = 1, lwd = 3)
lines(R, UPPER, lty = 2, lwd = 1)
segments(x0 = X, x1 = X, y0 = MIN - 100, y1 = MIN, lwd = 0.2)



dev.off()
########################################################################



########################################################################
### Figure 3: t-values for relative support 
########################################################################
K <- 60
T <- 10
MONTH <- 1:10
t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
    base <- village.data$base_3km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    
    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    
    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_0.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -5
COL.MAX <- 5
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################







########################################################################
### Figure 7: distribution of support estimates
########################################################################
pdf(file = "support_density.pdf", width = 6, height = 6)

A <- village.data$support.isaf - village.data$support.taliban
B <- village.data$support.isaf 
C <- village.data$support.taliban

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=x
cex.sub=x

min.X <- min(A, B, C)
max.X <- max(A, B, C)
max.Y <- 1


plot(density(A), xlab = "", ylab = "", main = "", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub, axes = F, xlim = c(min.X, max.X), ylim = c(0, max.Y), col = "black", lwd = 3)

par(new = T)

plot(density(B), xlab = "", ylab = "", main = "", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub, axes = F, xlim = c(min.X, max.X), ylim = c(0, max.Y), col = "red", lty = 2)

par(new = T)

plot(density(C), xlab = "Level of support", ylab = "Density", main = "", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub, xlim = c(min.X, max.X), ylim = c(0, max.Y), col = "blue", lty = 3)

legend("topleft", legend = c("ISAF - Taliban", "ISAF", "Taliban"), bty = "n", col = c("black", "red", "blue"), lwd = c(2, 1, 1), lty = c(1, 2, 3), cex = 1.2)

dev.off()
########################################################################
########################################################################



########################################################################
### Figure 8: Bivariate relationship between support and base
########################################################################
pdf(file = "support_base.pdf", width = 8, height = 8)

par(mfrow = c(2, 2))

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=x
cex.sub=x

x <- village.data$support.diff
y <- village.data$base_1km
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of bases"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 1km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)

x <- village.data$support.diff
y <- village.data$base_3km
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of bases"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 3km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)


x <- village.data$support.diff
y <- village.data$base_5km
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of bases"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 5km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)

x <- village.data$support.diff
y <- village.data$base_10km
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of bases"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 10km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)

dev.off()
############################################################



########################################################################
### Figure 9: Bivariate relationship between support and aid
########################################################################
pdf(file = "support_aid.pdf", width = 8, height = 8)

par(mfrow = c(2, 2))

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=x
cex.sub=x

x <- village.data$support.diff
y <- village.data$aid
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of aid projects"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 1km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)


x <- village.data$support.diff
y <- village.data$aid_3
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of aid projects"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 3km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)


x <- village.data$support.diff
y <- village.data$aid_5
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of aid projects"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 5km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)


x <- village.data$support.diff
y <- village.data$aid_10
XLAB <- "Relative level of ISAF support"
YLAB <- "Number of aid projects"
plot(x, y, xlab = XLAB, ylab = YLAB, main = "Within a 10km radius of villages", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
xx <- sort(x) 
lo <- loess(y ~ xx) 
lines(xx, predict(lo), col = "red", lwd = 2)

dev.off()

############################################################



########################################################################
### Figure 10: 15 km, 150 days (row data - simple bivariate correlation)
########################################################################
k <- 15
t <- 5
pre.A <- pre.ied.attack[[t]][, k]
pre.B <- pre.ied.found[[t]][, k]
pre.C <- pre.nonied[[t]][, k]
post.A <- post.ied.attack[[t]][, k]
post.B <- post.ied.found[[t]][, k]
post.C <- post.nonied[[t]][, k]
########################################################################
X<- village.data$support.isaf - village.data$support.taliban
########################################################################
pdf(file = "effect_15km_5month_correlation.pdf", width = 12, height = 4)

par(mfrow = c(1, 3))

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=2.2
cex.sub=x

MIN <- min(post.A, post.B, post.C)
MAX <- max(post.A, post.B, post.C)
XLAB = "Relative level of ISAF support"      
YLAB <- "Number of insurgent attacks"


MAIN <- "IED attacks"
plot(X, post.A, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN, MAX), col = 1, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(lowess(X, post.A), lwd = 2, col = "red")

MAIN <- "IED found"
plot(X, post.B, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN, MAX), col = 1, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(lowess(X, post.A), lwd = 2, col = "red")

MAIN <- "Non-IED attacks"
plot(X, post.C, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN, MAX), col = 1, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(lowess(X, post.A), lwd = 2, col = "red")



dev.off()
########################################################################


########################################################################
### Figure 11: standardized effects of relative support 
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

effect.A.sup <- effect.B.sup <- effect.C.sup <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
    base <- village.data$base_3km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    effect.A.sup[t, k] <- (coef(out)[2] * (max(X) - min(X))) / sd(Y)
        
    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    effect.B.sup[t, k] <- (coef(out)[2] * (max(X) - min(X))) / sd(Y)
    
    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    effect.C.sup[t, k] <- (coef(out)[2] * (max(X) - min(X))) / sd(Y)
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "std_effect_sup.png", width = 1800, height = 600)
effect.A <- effect.A.sup
effect.B <- effect.B.sup
effect.C <- effect.C.sup
COL.MIN <- -1.5
COL.MAX <- 1.5
BY <- 0.1

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)

Q <- effect.A
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- effect.B
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- effect.C
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################



########################################################################
### Figure 12: t-values for absolute level of ISAF support
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf 
    aid <- village.data$aid
    base <- village.data$base_3km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    
    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    
    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_isaf.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -5
COL.MAX <- 5
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################



########################################################################
### Figure 13: t-values for absolute level of taliban support
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.taliban 
    aid <- village.data$aid
    base <- village.data$base_3km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    
    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    
    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_taliban.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -5
COL.MAX <- 5
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################












########################################################################
### Figure 14: mahalanobis matching
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

verbose <- 1
t.value.A1 <- t.value.B1 <- t.value.C1 <- t.value.D1 <- t.value.E1 <- t.value.F1 <- matrix(NA, T, K) 
t.value.A2 <- t.value.B2 <- t.value.C2 <- t.value.D2 <- t.value.E2 <- t.value.F2 <- matrix(NA, T, K) 
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    N <- length(X)
    index <- 1:N
    aid <- village.data$aid
    base <- village.data$base_3km
########################################################################
    X.diff <- rep(NA, N)
    Y.diff <- rep(NA, N)
    Y.lag.diff <- rep(NA, N)
    base.diff <- rep(NA, N)
    aid.diff <- rep(NA, N)
    Y <- post.A
    Y.lag <- pre.A
    xxx <- cbind(Y.lag, base, aid)
    for (i in 1:N){
        index.i <- index[-i]
        xxx.i <- xxx[i, ]
        distance <- rep(0, N)
        for (j in 1:N){
            xxx.j <- xxx[j, ]
            distance[j] <- sqrt((xxx.i - xxx.j) %*% solve(cov(xxx)) %*% t(t(xxx.i) - xxx.j))
        }
        index.i <- index[-i]
        distance.i <- distance[-i]
        control <- index.i[distance.i == min(distance.i)]

        X.diff[i] <- X[i] - mean(X[control], na.rm = TRUE)
        Y.diff[i] <- Y[i] - mean(Y[control], na.rm = TRUE)
        Y.lag.diff[i] <- Y.lag[i] - mean(Y.lag[control], na.rm = TRUE)
        base.diff[i] <- base[i] - mean(base[control], na.rm = TRUE)
        aid.diff[i] <- aid[i] - mean(aid[control], na.rm = TRUE)
    }
    model <- lm(Y.diff ~ X.diff + Y.lag.diff + base.diff + aid.diff)
    beta <- coef(model)[2]
    se <- sqrt(vcovHC(model)[2, 2])
    A1 <- t.value.A1[t, k] <- beta/se
    
    X.diff <- rep(NA, N)
    Y.diff <- rep(NA, N)
    Y.lag.diff <- rep(NA, N)
    base.diff <- rep(NA, N)
    aid.diff <- rep(NA, N)
    Y <- post.B 
    Y.lag <- pre.B 
    xxx <- cbind(Y.lag, base, aid)
    for (i in 1:N){
        index.i <- index[-i]
        xxx.i <- xxx[i, ]
        distance <- rep(0, N)
        for (j in 1:N){
            xxx.j <- xxx[j, ]
            distance[j] <- sqrt((xxx.i - xxx.j) %*% solve(cov(xxx)) %*% t(t(xxx.i) - xxx.j))
        }
        index.i <- index[-i]
        distance.i <- distance[-i]
        control <- index.i[distance.i == min(distance.i)]

        X.diff[i] <- X[i] - mean(X[control], na.rm = TRUE)
        Y.diff[i] <- Y[i] - mean(Y[control], na.rm = TRUE)
        Y.lag.diff[i] <- Y.lag[i] - mean(Y.lag[control], na.rm = TRUE)
        base.diff[i] <- base[i] - mean(base[control], na.rm = TRUE)
        aid.diff[i] <- aid[i] - mean(aid[control], na.rm = TRUE)
    }
    model <- lm(Y.diff ~ X.diff + Y.lag.diff + base.diff + aid.diff)
    beta <- coef(model)[2]
    se <- sqrt(vcovHC(model)[2, 2])
    B1 <- t.value.B1[t, k] <- beta/se

    X.diff <- rep(NA, N)
    Y.diff <- rep(NA, N)
    Y.lag.diff <- rep(NA, N)
    base.diff <- rep(NA, N)
    aid.diff <- rep(NA, N)
    Y <- post.C 
    Y.lag <- pre.C 
    xxx <- cbind(Y.lag, base, aid)
    for (i in 1:N){
        index.i <- index[-i]
        xxx.i <- xxx[i, ]
        distance <- rep(0, N)
        for (j in 1:N){
            xxx.j <- xxx[j, ]
            distance[j] <- sqrt((xxx.i - xxx.j) %*% solve(cov(xxx)) %*% t(t(xxx.i) - xxx.j))
        }
        index.i <- index[-i]
        distance.i <- distance[-i]
        control <- index.i[distance.i == min(distance.i)]

        X.diff[i] <- X[i] - mean(X[control], na.rm = TRUE)
        Y.diff[i] <- Y[i] - mean(Y[control], na.rm = TRUE)
        Y.lag.diff[i] <- Y.lag[i] - mean(Y.lag[control], na.rm = TRUE)
        base.diff[i] <- base[i] - mean(base[control], na.rm = TRUE)
        aid.diff[i] <- aid[i] - mean(aid[control], na.rm = TRUE)
    }
    model <- lm(Y.diff ~ X.diff + Y.lag.diff + base.diff + aid.diff)
    beta <- coef(model)[2]
    se <- sqrt(vcovHC(model)[2, 2])
    C1 <- t.value.C1[t, k] <- beta/se

    if (verbose!=0&k%%verbose == 0){
      cat("t", t, "k", k, "A1", A1, "B1", B1, "C1", C1, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_match_0.png", width = 1800, height = 600)


ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -4.5
COL.MAX <- 4.5
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)


dev.off()
##############################################################



########################################################################
### Figure 15: t-value of relative support: Base = 1, 5, 10km
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)
t.value.A4 <- t.value.B4 <- t.value.C4 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
########################################################################
    base <- village.data$base_1km
    #base <- village.data$base_3km
    #base <- village.data$base_5km
    #base <- village.data$base_10km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.A1[t, k] <- coef.out / se.out

    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.B1[t, k] <- coef.out / se.out

    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.C1[t, k] <- coef.out / se.out
########################################################################
    #base <- village.data$base_1km
    #base <- village.data$base_3km
    base <- village.data$base_5km
    #base <- village.data$base_10km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.A2[t, k] <- coef.out / se.out

    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.B2[t, k] <- coef.out / se.out

    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.C2[t, k] <- coef.out / se.out
########################################################################
    #base <- village.data$base_1km
    #base <- village.data$base_3km
    #base <- village.data$base_5km
    base <- village.data$base_10km
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.A3[t, k] <- coef.out / se.out

    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.B3[t, k] <- coef.out / se.out

    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.C3[t, k] <- coef.out / se.out
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_b.png", width = 2200, height = 2200) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -5
COL.MAX <- 5
BY <- 0.5




Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- t.value.A2
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B2
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C2
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P6 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.A3
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P7 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B3
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P8 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C3
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P9 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )



ad1 <- 0.04
ad2 <- -0.02
ad3 <- 0.0
print(P1, position = c(0, 0.69 + ad2, 0.315, 0.97 -ad1), more = TRUE)
print(P2, position = c(0.315, 0.69  + ad2, 0.63, 0.97 -ad1), more = TRUE)
print(P3, position = c(0.63, 0.69  + ad2, 1 , 0.97 -ad1), more = TRUE)
 
print(P4, position = c(0, 0.35 + ad2, 0.315, 0.64 -ad1), more = TRUE)
print(P5, position = c(0.315, 0.35  + ad2, 0.63 , 0.64 -ad1), more = TRUE)
print(P6, position = c(0.63 , 0.35  + ad2, 1, 0.64 -ad1), more = TRUE)
 
print(P7, position = c(0, 0.02 + ad2, 0.315, 0.30 -ad1), more = TRUE)
print(P8, position = c(0.315, 0.02  + ad2 , 0.63, 0.30 -ad1), more = TRUE)
print(P9, position = c(0.63 , 0.02  + ad2 , 1, 0.30 -ad1))

panel.text(x = 1050, y = 95, labels = "Bases within 1km", cex = 4, font = 1)

panel.text(x = 1050, y = 825, labels = "Bases within 5km", cex = 4, font = 1)

panel.text(x = 1050, y = 1575, labels = "Bases within 10km", cex = 4, font = 1)



dev.off()
################################################################################



########################################################################
### Figure 16: t-value for relative support: Aid = 3, 5, 10km
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)
t.value.A4 <- t.value.B4 <- t.value.C4 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    base <- village.data$base_3km
########################################################################
    aid <- village.data$aid_3
    #aid <- village.data$aid_5
    #aid <- village.data$aid_10
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.A1[t, k] <- coef.out / se.out

    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.B1[t, k] <- coef.out / se.out

    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.C1[t, k] <- coef.out / se.out
########################################################################
    #aid <- village.data$aid_3
    aid <- village.data$aid_5
    #aid <- village.data$aid_10
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.A2[t, k] <- coef.out / se.out

    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.B2[t, k] <- coef.out / se.out

    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.C2[t, k] <- coef.out / se.out
########################################################################
    #aid <- village.data$aid_3
    #aid <- village.data$aid_5
    aid <- village.data$aid_10
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.A3[t, k] <- coef.out / se.out

    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.B3[t, k] <- coef.out / se.out

    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    t.value.C3[t, k] <- coef.out / se.out
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_a.png", width = 2200, height = 2200) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -5
COL.MAX <- 5
BY <- 0.5




Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- t.value.A2
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B2
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C2
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P6 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.A3
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P7 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B3
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P8 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C3
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P9 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )




ad1 <- 0.04
ad2 <- -0.02
ad3 <- 0.0
print(P1, position = c(0, 0.69 + ad2, 0.315, 0.97 -ad1), more = TRUE)
print(P2, position = c(0.315, 0.69  + ad2, 0.63, 0.97 -ad1), more = TRUE)
print(P3, position = c(0.63, 0.69  + ad2, 1 , 0.97 -ad1), more = TRUE)
 
print(P4, position = c(0, 0.35 + ad2, 0.315, 0.64 -ad1), more = TRUE)
print(P5, position = c(0.315, 0.35  + ad2, 0.63 , 0.64 -ad1), more = TRUE)
print(P6, position = c(0.63 , 0.35  + ad2, 1, 0.64 -ad1), more = TRUE)
 
print(P7, position = c(0, 0.02 + ad2, 0.315, 0.30 -ad1), more = TRUE)
print(P8, position = c(0.315, 0.02  + ad2 , 0.63, 0.30 -ad1), more = TRUE)
print(P9, position = c(0.63 , 0.02  + ad2 , 1, 0.30 -ad1))


panel.text(x = 1050, y = 95, labels = "Aid projects within 3km", cex = 4, font = 1)

panel.text(x = 1050, y = 825, labels = "Aid projects within 5km", cex = 4, font = 1)

panel.text(x = 1050, y = 1575, labels = "Aid projects within 10km", cex = 4, font = 1)



dev.off()
################################################################################




########################################################################
### Figure 17: t-values for ISAF control 
########################################################################
K <- 60
T <- 10
MONTH <- 1:10
t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
    base <- village.data$base_3km
    control <- village.data$control
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    ##out <- lm(Y ~ X + Y.lag + base + aid + control)
    out <- lm(Y ~ control + Y.lag + base + aid + X)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    
    Y <- post.B 
    Y.lag <- pre.B 
    ##out <- lm(Y ~ X + Y.lag + base + aid + control)
    out <- lm(Y ~ control + Y.lag + base + aid + X)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    
    Y <- post.C 
    Y.lag <- pre.C 
    ##out <- lm(Y ~ X + Y.lag + base + aid + control)
    out <- lm(Y ~ control + Y.lag + base + aid + X)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_isaf_control_effect_of_control.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -8
COL.MAX <- 8
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################




########################################################################
### Figure 18: t-values for relative support with ISAF control as a covariate
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
    base <- village.data$base_3km
    control <- village.data$control
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid + control)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    
    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid + control)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    
    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid + control)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_isaf_control_effect_of_support.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -6
COL.MAX <- 6
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################




########################################################################
### Figure 19: distribution of ideal point estimates 
########################################################################
pdf(file = "ideal_point_density.pdf", width = 6, height = 6)

A <- village.data$ideal.point

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=x
cex.sub=x

plot(density(A), xlab = "Ideal point for policy reform", ylab = "Density", main = "", cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub, col = "black", lwd = 3)


dev.off()
########################################################################
########################################################################



########################################################################
### Figure 20: t-values for relative support with ideal point estimates as a covariate
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
    base <- village.data$base_3km
    ideal <- village.data$ideal.point
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out <- lm(Y ~ X + Y.lag + base + aid + ideal)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    
    Y <- post.B 
    Y.lag <- pre.B 
    out <- lm(Y ~ X + Y.lag + base + aid + ideal)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    
    Y <- post.C 
    Y.lag <- pre.C 
    out <- lm(Y ~ X + Y.lag + base + aid + ideal)
    coef.out <- coef(out)[2]
    se.out <- sqrt(vcovHC(out)[2, 2])
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_with_ideal_point_effect_of_support.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -5
COL.MAX <- 5
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################




########################################################################
### Figure 21: t-values for relative support with district-level random effects
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- matrix(NA, T, K)
t.value.A3 <- t.value.B3 <- t.value.C3 <- matrix(NA, T, K)
PA <- PB <- PC <- matrix(NA, T, K)

verbose <- 1
for (t in 1:T){
  for (k in 1:K){
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
########################################################################
    X<- village.data$support.isaf - village.data$support.taliban
    aid <- village.data$aid
    base <- village.data$base_3km
    dist <- village.data$d398distid
    prov <- village.data$prov
########################################################################
    ### robust se ###
    Y <- post.A
    Y.lag <- pre.A
    out.lm <- lm(Y ~ X + Y.lag + base + aid)
    out.lmer <- lmer(Y ~ X + Y.lag + base + aid + (1|dist), REML = TRUE)
    out <- out.lmer
    coef.out <- fixef(out)[2]
    se.out <- sqrt(vcov(out)[2, 2])    
    A1 <- t.value.A1[t, k] <- coef.out / se.out
    A2 <- t.value.A2[t, k] <- coef.out
    A3 <- t.value.A3[t, k] <- se.out
    PA[t, k] <- anova(out.lmer, out.lm)[[8]][2]
    
    Y <- post.B 
    Y.lag <- pre.B 
    out.lm <- lm(Y ~ X + Y.lag + base + aid)
    out.lmer <- lmer(Y ~ X + Y.lag + base + aid + (1|dist), REML = TRUE)
    out <- out.lmer
    coef.out <- fixef(out)[2]
    se.out <- sqrt(vcov(out)[2, 2])   
    B1 <- t.value.B1[t, k] <- coef.out / se.out
    B2 <- t.value.B2[t, k] <- coef.out
    B3 <- t.value.B3[t, k] <- se.out
    PB[t, k] <- anova(out.lmer, out.lm)[[8]][2]
    
    Y <- post.C 
    Y.lag <- pre.C 
    out.lm <- lm(Y ~ X + Y.lag + base + aid)
    out.lmer <- lmer(Y ~ X + Y.lag + base + aid + (1|dist), REML = TRUE)
    out <- out.lmer
    coef.out <- fixef(out)[2]
    se.out <- sqrt(vcov(out)[2, 2])   
    C1 <- t.value.C1[t, k] <- coef.out / se.out
    C2 <- t.value.C2[t, k] <- coef.out
    C3 <- t.value.C3[t, k] <- se.out
    PC[t, k] <- anova(out.lmer, out.lm)[[8]][2]
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
    
  }
}
########################################################################
png(file = "t_value_dist_random_effect.png", width = 1800, height = 600) 

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
##COL.MIN <- -5
##COL.MAX <- 5
COL.MIN <- -3.2
COL.MAX <- 3.2
BY <- 0.5


Q <- t.value.A1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.B1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

BY <- 0.2
Q <- t.value.C1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P3, position = c(0.63, 0, 1, 1), more = TRUE)

dev.off()
################################################################################





########################################################################
### Figure 22: correlation between districts and covariates 
########################################################################
pdf("correlation_district_covariates.pdf", height = 6, width = 10)
########################################################################
k <- 15
t <- 5
pre.A <- pre.ied.attack[[t]][, k]
pre.B <- pre.ied.found[[t]][, k]
pre.C <- pre.nonied[[t]][, k]
########################################################################
support <- village.data$support.isaf - village.data$support.taliban
aid <- village.data$aid
base <- village.data$base_3km
dist <- village.data$d398distid
prov <- village.data$prov
dist.prov <- data.frame(dist, prov)
########################################################################
### Logar(5); Kunar(10); Helmand(23); Uruzgan(26); Khost(32)
########################################################################

par(mfrow = c(2, 3))
for (h in 1:6){
    if(h == 1){
        MAIN <- "Relative ISAF support"
        x <- tapply(support, dist, mean)
        ADJ <- 0.5
    }
    if(h == 2){
        MAIN <- "Number of ISAF bases"
        x <- tapply(base, dist, mean)
        ADJ <- 0.5
    }
    if(h == 3){
        MAIN <- "Number of aid projects"
        x <- tapply(aid, dist, mean)
        ADJ <- 0.5
    }
    if(h == 4){
        MAIN <- "Number of past IED attacks \n (15km, 150 days)"
        x <- tapply(pre.A, dist, mean)
        ADJ <- 20
    }
    if(h == 5){
        MAIN <- "Number of past IED found \n (15km, 150 days)"
        x <- tapply(pre.B, dist, mean)
        ADJ <- 30
    }
    if(h == 6){
        MAIN <- "Number of past non-IED attacks \n (15km, 150 days)"
        x <- tapply(pre.C, dist, mean)
        ADJ <- 150
    }
    name <- as.numeric(rownames(x))
    MIN <- min(x)
    MAX <- max(x)
    plot(x, cex = 0, axes = F, xlab = "", ylab = "", xlim = c(1, length(x)), ylim = c(MIN, MAX + ADJ))
    par(new = T)
    for (i in 1:length(x)){
        prov.i <- unique(dist.prov$prov[dist.prov$dist == name[i]])
        if(prov.i == 5){
            COL <- "black"
            PCH <- 1
        }
        if(prov.i == 10){
            COL <- "red"
            PCH <- 2
        }
        if(prov.i == 23){
            COL <- "blue"
            PCH <- 3
        }
        if(prov.i == 26){
            COL <- "orange"
            PCH <- 4
        }
        if(prov.i == 32){
            COL <- "purple"
            PCH <- 5
        }
        
        plot(i, x[i], xlab = "", ylab = "", xlim = c(1, length(x)), ylim = c(MIN, MAX + ADJ), pch = PCH, col = COL, axes = F)
        par(new = T)
    }
    plot(x, cex = 0, xlab = "District IDs", ylab = "", xlim = c(1, length(x)), ylim = c(MIN, MAX + ADJ), main = MAIN)
    
    text(x = c(2, 5.5, 10, 15, 19), y = rep(MAX + ADJ, 5), labels = c("Logar", "Kunar", "Helmand", "Uruzgan", "Khost"), col = c("black", "red", "blue", "orange", "purple"), cex = rep(1.3, 5))
}

dev.off()
########################################################################
########################################################################





########################################################################
### Figure 27: Negative Binomial Model: 15 km, 150 days
########################################################################
k <- 15
t <- 5

pre.A <- pre.ied.attack.civilian[[t]][, k]
pre.C <- pre.nonied.civilian[[t]][, k]
post.A <- post.ied.attack.civilian[[t]][, k]
post.C <- post.nonied.civilian[[t]][, k]
pre.E <- pre.A + pre.C
post.E <- post.A + post.C
########################################################################
J <- 100
mean1 <- mean2 <- mean3 <- mean4 <- mean5 <- mean6 <- mean7 <- mean8 <- mean9 <- mean10 <- mean11 <- mean12 <- rep(NA, J)
lower1 <- lower2 <- lower3 <- lower4 <- lower5 <- lower6 <- lower7 <- lower8 <- lower9 <- lower10 <- lower11 <- lower12 <- rep(NA, J)
upper1 <- upper2 <- upper3 <- upper4 <- upper5 <- upper6 <- upper7 <- upper8 <- upper9 <- upper10 <- upper11 <- upper12 <- rep(NA, J)
########################################################################
X<- village.data$support.isaf - village.data$support.taliban
aid <- village.data$aid
base <- village.data$base_3km
MIN <- min(X, na.rm = TRUE)
MAX <- max(X, na.rm = TRUE)
R <- seq(MIN, MAX, length.out = J)
########################################################################
### robust se ###
Y <- post.A
Y.lag <- pre.A
out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500))
coef.out <- coef(out)
vcov.out <- vcovHC(out)
set.seed(1000)
beta <- mvrnorm(1000, mu = coef.out, Sigma = vcov.out)
mean.Y.lag <- mean(Y.lag, na.rm = TRUE)
mean.base <- mean(base, na.rm = TRUE)
median.aid <- median(aid, na.rm = TRUE)
Z <- cbind(1, R, mean.Y.lag, mean.base, median.aid)  
mu <- exp(beta %*% t(Z))
for (j in 1:J){
  mu.j <- mu[, j]
  mean1[j] <- mean(mu.j)
  upper1[j] <- quantile(mu.j, prob = 0.975)
  lower1[j] <- quantile(mu.j, prob = 0.025)    
}

Y <- post.C 
Y.lag <- pre.C 
out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500))
coef.out <- coef(out)
vcov.out <- vcovHC(out)
set.seed(1000)
beta <- mvrnorm(1000, mu = coef.out, Sigma = vcov.out)
mean.Y.lag <- mean(Y.lag, na.rm = TRUE)
mean.base <- mean(base, na.rm = TRUE)
median.aid <- median(aid, na.rm = TRUE)
Z <- cbind(1, R, mean.Y.lag, mean.base, median.aid)   
mu <- exp(beta %*% t(Z))
for (j in 1:J){
  mu.j <- mu[, j]
  mean3[j] <- mean(mu.j)
  upper3[j] <- quantile(mu.j, prob = 0.975)
  lower3[j] <- quantile(mu.j, prob = 0.025)    
}

Y <- post.E
Y.lag <- pre.E
out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500))
coef.out <- coef(out)
vcov.out <- vcovHC(out)
set.seed(1000)
beta <- mvrnorm(1000, mu = coef.out, Sigma = vcov.out)
mean.Y.lag <- mean(Y.lag, na.rm = TRUE)
mean.base <- mean(base, na.rm = TRUE)
median.aid <- median(aid, na.rm = TRUE)
Z <- cbind(1, R, mean.Y.lag, mean.base, median.aid)   
mu <- exp(beta %*% t(Z))
for (j in 1:J){
  mu.j <- mu[, j]
  mean5[j] <- mean(mu.j)
  upper5[j] <- quantile(mu.j, prob = 0.975)
  lower5[j] <- quantile(mu.j, prob = 0.025)    
}

########################################################################
### computing change in attacks
########################################################################
RR <- c(-2.5, 0.5)
Z <- cbind(1, RR, mean.Y.lag, mean.base, median.aid)   
mu <- beta %*% t(Z)
diff <- mu[,2] - mu[,1]
quantile(diff, prob = 0.025)
mean(diff)
quantile(diff, prob = 0.975)
########################################################################
pdf(file = "effect_15km_5month_civilian_nb.pdf", width = 12, height = 4)

par(mfrow = c(1,3))

#################################
X<- village.data$support.isaf - village.data$support.taliban
XLAB = "Relative level of ISAF support"      
YLAB <- "Number of insurgent attacks"
ADJ.FACTOR <- 50
#################################


MIN <- min(lower1, lower3, lower5)
MAX <- max(upper1, upper3, upper5)

x <- 1.5
cex.lab=x
cex.axis=x
cex.main=2.2
cex.sub=x

LOWER <- lower1
MEAN <- mean1
UPPER <- upper1
MAIN <- "IED attacks"
ADJ <- abs(MAX - MIN)/ADJ.FACTOR
plot(R, MEAN, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN - ADJ, MAX), col = 0, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(R, LOWER, lty = 2, lwd = 1)
lines(R, MEAN, lty = 1, lwd = 3)
lines(R, UPPER, lty = 2, lwd = 1)
segments(x0 = X, x1 = X, y0 = MIN - 100, y1 = MIN, lwd = 0.2)

LOWER <- lower3
MEAN <- mean3
UPPER <- upper3
MAIN <- "Non-IED attacks"
ADJ <- abs(MAX - MIN)/ADJ.FACTOR
plot(R, MEAN, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN - ADJ, MAX), col = 0, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(R, LOWER, lty = 2, lwd = 1)
lines(R, MEAN, lty = 1, lwd = 3)
lines(R, UPPER, lty = 2, lwd = 1)
segments(x0 = X, x1 = X, y0 = MIN - 100, y1 = MIN, lwd = 0.2)

LOWER <- lower5
MEAN <- mean5
UPPER <- upper5
MAIN <- "Both"
ADJ <- abs(MAX - MIN)/ADJ.FACTOR
plot(R, MEAN, xlab = XLAB, ylab = YLAB, main = MAIN, ylim = c(MIN - ADJ, MAX), col = 0, cex.lab = cex.lab, cex.axis = cex.axis, cex.main = cex.main, cex.sub = cex.sub)
lines(R, LOWER, lty = 2, lwd = 1)
lines(R, MEAN, lty = 1, lwd = 3)
lines(R, UPPER, lty = 2, lwd = 1)
segments(x0 = X, x1 = X, y0 = MIN - 100, y1 = MIN, lwd = 0.2)



dev.off()
########################################################################


########################################################################
### Figure 28: t-values (Negative Binomial Model)
########################################################################
K <- 60
T <- 10
MONTH <- 1:10

t.value.A1 <- t.value.B1 <- t.value.C1 <- t.value.D1 <- t.value.E1  <- t.value.F1 <- matrix(NA, T, K)
t.value.A2 <- t.value.B2 <- t.value.C2 <- t.value.D2 <- t.value.E2  <- t.value.F2 <- matrix(NA, T, K)

verbose <- 1

T.min <- 1
K.min <- 1

for (t in T.min:T){
    for (k in K.min:K){
        pre.D <- pre.ied.attack.civilian[[t]][, k]
        pre.E <- pre.nonied.civilian[[t]][, k]
        post.D <- post.ied.attack.civilian[[t]][, k]
        post.E <- post.nonied.civilian[[t]][, k]
        pre.F <- pre.D + pre.E
        post.F <- post.D + post.E
########################################################################
        X<- village.data$support.isaf - village.data$support.taliban
        aid <- village.data$aid
        base <- village.data$base_3km
########################################################################
        #options(warn = 2)
        tt <- t
########################################################################
### robust se ###
        out <- NULL
        D1 <- NULL
        Y <- post.D 
        Y.lag <- pre.D 
        try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)), silent = TRUE)
        if(!is.null(out)){
            coef.out <- coef(out)[2]
            se.out <- sqrt(vcovHC(out)[2, 2])
            D1 <- t.value.D1[tt, k] <- coef.out / se.out
        }

        out <- NULL
        E1 <- NULL
        Y <- post.E 
        Y.lag <- pre.E 
        try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)), silent = TRUE)
        if(!is.null(out)){
            coef.out <- coef(out)[2]
            se.out <- sqrt(vcovHC(out)[2, 2])
            E1 <- t.value.E1[tt, k] <- coef.out / se.out
        }

        out <- NULL
        F1 <- NULL
        Y <- post.F 
        Y.lag <- pre.F 
        try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)), silent = TRUE)
        if(!is.null(out)){
            coef.out <- coef(out)[2]
            se.out <- sqrt(vcovHC(out)[2, 2])
            F1 <- t.value.F1[tt, k] <- coef.out / se.out
        }
        
        if (verbose!=0&t%%verbose == 0){
            cat("t", t, "k", k, "D1", D1, "E1", E1, "F1", F1, "\n") 
        }   
    }
}
########################################################################
png(file = "t_value_civilian_nb.png", width = 1800, height = 600)

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -8
COL.MAX <- 8
BY <- 0.5

Q <- t.value.D1
MAIN <- list("IED attacks", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.E1
MAIN <- list("Non-IED attacks", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                #colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- t.value.F1
MAIN <- list("Both", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )




print(P1, position = c(0, 0, 0.315, 1), more = TRUE)
print(P3, position = c(0.315, 0, 0.63, 1), more = TRUE)
print(P5, position = c(0.63, 0, 1, 1), more = TRUE)


dev.off()
################################################################################

